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Abstract 

Molecular dynamics simulations are used to investigate the atomic mobil- 
ity and diffusivity of a generalized Frenkel-Kontorova model which takes into 
account anharmonic (exponential) interaction of atoms subjected to a three- 
dimensional substrate potential periodic in two dimensions and nonconvex 
(Morse) in the third dimension. The numerical results are explained by a 
phenomenological theory which treats a system of strongly interacting atoms 
as a system of weakly interacting quasiparticles (kinks). Model parameters 
are chosen close to those for K-W(112) adsorption system. 
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I. INTRODUCTION 



Experimental studies of transport coefficients in systems of strongly interacting atoms 
adsorbed on a crystalline surface show a very rich and complicated behavior, especially 
as functions of the atomic concentration. The variation of the diffusion coefficient versus 
coverage is particularly important for adsorbed layers where the concentration may be varied 
in wide limits from zero (diffusion of isolated adsorbed atoms) to very high values (for 
example, in some adsystems the interatomic distance in a monolayer of adatoms is lower 
than that in the corresponding massive crystal)0. The theoretical study of mass and charge 
transport in such systems is a very difficult problem; however it was studied for various 
kinds of interactions by Gomer et all using Monte Carlo simulations. At high temperatures, 
transport coefficients can be found with a perturbation technique starting from the case 
of noninteracting atomsi. At low temperatures, the case of interacting atoms has been 
studied by a numerical calculation of the transport properties of the standard Frenkel- 
Kontorova (FK) model, which describes a chain of harmonically interacting atoms subjected 



of a system of strongly interacting atoms in a more general one- dimensional model has 
been approximately treated with a phenomenological approach which introduces weakly 
interacting quasi-particle^l. This method provides analytical estimates for the transport 
coefficients, but it requires many approximations. In particular the properties of the quasi- 
particles involved in the theory are deduced from results of the standard FK model which 
provides only a simplified picture. Therefore it was necessary to check the validity of the 
theoretical approach by numerical simulations of a model which is sufficiently complicated to 
provide a reasonable description of a real system. This is the aim of this paper which studies 
a two-dimensional generalized FK model and also discusses some experimental results in the 
same perspective. 

The original FK model was introduced to analyze the dynamics of dislocations in crystal^l 
by considering a chain of interacting particles subjected to a periodic substrate (on-site) 




Recently the low-temperature behavior 
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potential. It can describe, for example, a closely-packed row of atoms in a crystaEi, a chain 
of atoms adsorbed on stepped or furrowed crystal surfaced, a chain of ions in a "channel" 
of a quasi-one-dimensional conductoiS, hydrogen atoms in hydrogen-bonded systems^, etc. 
In all the cases mentioned above, the chain of interacting particles is an intrinsic part of the 
whole physical system under consideration. The role of the remainder of the system is played 
by an external substrate potential and a thermal bath. Although it is still oversimplified, the 
generalized FK model that we consider here provides a rather complete description of a layer 
of atoms adsorbed on a two-dimensional crystal surface. It includes realistic (exponential) 
interactions of particles instead of the harmonic springs of the standard FK model, and the 
substrate potential is three-dimensional. It is periodic in the two dimensions parallel to the 
surface and has a Morse shape in the third direction, orthogonal to the surface. 

The transport properties of the system are described by two coefficients, the mobility 
B and the chemical diffusivity Dc- The mobility defines the response of the system to an 
infinitesimal d.c. force F, 

J = pBF, (1) 

where J is the atomic fiux caused by the force and p is the average atomic concentration. On 
the other hand, the chemical diffusion coefficient Dc connects the flux J{x, t) in a nonequi- 
librium state to the gradient of the atomic concentration when p{x, t) slightly deviates from 
its equilibrium value. According to Fick's law 

{{J{x,t)))^-Dc-^{{p{x,t))), (2) 

where ((. . .)) stands for the averaging over macroscopic distances x ^ a^, and is the 
average interatomic distance. These two coefficients are coupled through the relation 

^ (3) 

X 

where /c^ is Boltzmann's constant, T the temperature and x the dimensionless susceptibility 
of the system. 



The predictions of the phenomenological approachQ can be summarized as follows. The 
mass transport is caused by kinks which describe localized compressions or expansions of 
the chain and therefore the mobility B can be expected to be proportional to the kink 
concentration. The kinks have two different origins, "geometrical" and thermal. We call 
geometrical kinks, the kinks which result from the value of the coverage 6 = N/M, where 

is the number of atoms and M the number of wells of the substrate potential on a given 
length. For 6 = 1/q, with integer g (g = 1, 2, . . .), the system has a trivial ground state with 
one atom at the bottom of the substrate wells every q^^ well. When 9 deviates slightly from 
such a value, the difference is accommodated by the system by forming localized discom- 
mensurations which are the geometrical kinks (called also "trivial kinks" in the notation of 
Ref. P). As the kink density increases when 6 deviates from 1/q, the theory predicts that 
B{6) should exhibit local minima for any trivial ground state (GS) of the system, such as 
6 = 1,6 = 1/2, 6 = 1/3, etc. When 6 = p/q is a rational number with a larger numerator, 
such as 6* = 2/3, the density of geometrical kinks becomes very large and one could expect 
to get a high mobility B. The picture is however more complicated because, due to their 
high density, the geometrical kinks interact strongly and, when temperature is sufficiently 
low with respect to their interaction energy, they tend to form a regular lattice which is 
weakly pinned, giving a low mobility for any rational 6. A slight deviation from 6 = p/q 
appears as discommensurations in the kink lattice i.e. "kinks in a kink lattice", which are 
called superkinks in Ref. ^. These topological excitations of the kink lattice contribute to 
mass transport exactly as the trivial kink do, so that the mobility is expected to exhibit 
local minima for 6 = p/q such as 6 = 2/3. In the limit T — the function B{6) should 
therefore have minima at any rational 6. When temperature increases, the secondary min- 
ima disappear because the kink lattice "melts" and moreover thermal fluctuations create 
kink-antikink pairs which are thermally activated. Consequently at high enough temper- 
ature the mobility is expected to exhibit broad maxima between the primary minima at 
6 = 1/q. Such a behavior has been observed in one dimensional mo dels&H. 

The behavior of the diffusion coefficient Dc is simpler than the variation of B{6) as pre- 
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dieted in Ref. |T5|. Aeeording to Fiek's law (0), is the proportionality coefficient between 
the (infinitesimal) gradient of the atomic concentration and the fiux of atoms caused by this 
gradient. However, a gradient of atomic concentration automatically produces a correspond- 
ing gradient of kink concentration. In the standard FK model, where the elastic constant 
does not depend on 9 and where the parameters of kinks and antikinks are the same, Dc{6) is 
the ratio of two quantities which vary similarly so that it should be approximately constant 
and coincide with the kink (or antikink) diffusion coefficient. In the generalized FK model 
the situation is different because the anharmonicity of the interatomic interaction destroys 
the kink-antikink symmetryEl. The effective interatomic forces for a kink, which corresponds 
to a local contraction of the chain, exceed those for an antikink which is associated to region 
of a local extension. Thus, in comparison with an antikink, a kink is characterized by a larger 
value of the rest energy and by lower values of effective mass and activation energy for its 
motionEl. When the coverage passes through a commensurate value ^o, the geometrical-kink 
density vanishes; for 6 < Oq the system has geometrical antikinks while for 9 > Oq the system 
has geometrical kinks. Therefore, when the coverage 6 increases through a commensurate 
value 6q, the activation energy for the chemical diffusion should jump to a smaller value. 
Simultaneously the value of Dc should rise sharply when the coverage 6 exceeds the value 
6o that characterizes a "well-defined" commensurate structure and one could expect Dc{9) 
to exhibit the shape of a Devil's staircase. The abrupt (jump-like) increase of Dc{6) will 
only exist in the T ^ limit and, for any T 7^ 0, these jumps will be smoothed owing to 
corrections from thermally excited kink-antikink pairs. 

In present paper we check these predictions by molecular dynamics investigations of 
the low-temperature mobility and diffusivity of a generalized FK model in one and two 
dimensions. In Sec. ^ we describe the model and define its parameters. Kink parameters 
are calculated in Sec. |T|. Simulation results for the mobility are presented in Sec. 0, 
and those for the chemical diffusivity are described in Sec. [V|. Sec. |V| discusses known 
experimental results in the framework of these studies and Sec. |Vli] concludes the paper. 
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II. THE MODEL 



As for the standard FK model, we consider the dynamics of atoms adsorbed on a periodic 
substrate. The displacement of each atom is characterized by three variables: x and y 
describe its motion parallel to the surface, while z describes its deviation orthogonal to the 
substrate. For the substrate potential, we take the function 

Vsub{x, y, z) = [Vpr{x; a^x, Ssx, s^) + Vpr{y] asy, Ssy, Sy)]e~'^''' + Vz{z). (4) 

To model the substrate potential along the surface, we use a deformable periodic potential 
which can be adjusted to describe an actual crystal field0, 

V (^.n r X _ (1 + Sxf[^ - cos(27ra:/a,x)] 

Via;, asx, e^x, s,) - ^ ^ + ^2 _ cos(27rx/a,,.) ' ^ ^ 

Thus, Esx corresponds to the activation energy for diffusion of a single atom along the x 
direction, asx to the lattice constant and the parameter Sx {\sx\ < 1) controls the shape of 
the substrate potential. The frequency of a single-atom vibration along the x direction 
is connected to the shape parameter by the relationship = uJo{l + Sx)/{1 — Sx), where 
uiQ = {esx/'^'raY^'^{2n / asx) and m is the atomic mass. The potential Vpriy'idsyi^sy, Sy) has 
the same form. 

The potential perpendicular to the surface is modeled by the Morse function 

V,{z)=ed(e-^-^ (6) 

which tends to the adsorption energy Sd when z goes to infinity. The anharmonicity param- 
eter 7 is related to the frequency Uz of a single-atom vibration in the normal direction by 
the relation ojI = 2j'^ed/m. 

Finally, the exponential factor after the square brackets of the right-hand side of Eq. @ 
takes into account the decrease of the influence of the surface corrugation as the atoms move 
away from the surface, so that Vsubix, y, z) Sd when z 00. 

For the interaction between the atoms we take the exponential repulsion 
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V,nt{r) = VoeM-M, (7) 

where Vq is the amphtude and (3q^ determines the typical range of the interaction. This 
potential is adapted to describe rather high coverages such that the atoms interact through 
the repulsive branch of the interatomic potential. In numerical simulations, we can only 
include the interaction of a given adatom with a finite number of neighbors. Therefore, 
we use the standard approach of MD simulations and introduce a cutoff distance r*. We 
account only for the interactions between the atoms separated by distances lower than r* 
and to reduce errors caused by this procedure, the interaction potential is truncated as 

VUr) = V,nt{r) - Vintin - VUr*)ir - r*), (8) 

so that the interaction potential and force vanish at the cutoff distance, Vint{r*) = V/^^{r*) = 
(tilde will be omitted in what follows). Besides, because we are using the repulsive in- 
teratomic interaction, we have to fix the atomic concentration. It is imposed by periodic 
boundary conditions in x and y. We place atoms in the fixed area Lx x Ly, Lx = M^agx, 
Ly = Myttsy, so that the dimensionless atomic concentration is equal to 6* = N/M, where 

M = MxMy. 

To model the energy exchange of the atoms with a thermal bath, we use the Langevin 
equations for atomic coordinates Xi 



d 

mxi + mrjXi + - — 
dxi 



Vsubixi, yi, Zi) + Vint{\ri - fj\) 



F'^-''> + 6Ft\t), (9) 



and similar equations for y^ and Zi. Here, rj corresponds to the rate of the energy exchange 
with the substrate, F = {F, 0, 0} to the dc driving force, and 6F is a Gaussian random force 
with correlation function 

5Ff = 2r,mkBT5^p5,,5{t - t'). (10) 

We use a dimensional system of units adapted to the scales of the problem. Distance is 
measured in Angstroms, energy and temperature in eV. The mass of an adatom is chosen 
as our mass unit (m = 1). This imposes a time scale. We measure time in units of the 



characteristic time interval to = '^'^/^x- In the remainder of the paper, the measures of 
other dimensional physical quantities will be omitted, but they are all expressed in terms of 
the above-defined units. 

In order to be close to real physical systems, let us take the adsystem K-W(112) as 
an example to define the model parameters: the W(112) surface is characterized by a 
strong anisotropy of the atomic relief because it has close packed rows of substrate atoms 
separated by furrows of atomic depth. Namely, in the simulation, we put asx = 2.74 A and 
asy = 4.47 A which are respectively the distances between the neighboring wells along and 
across the furrows on the W(112) surface, and Sgx = 0.46 eV and Sgy = 0.76 eV for the 
corresponding barriers (these values were taken from Ref. |18]). To model the shape of the 



substrate potential, we have to define the parameters and Sy. They can be estimated 
to be within the range [0.2,0.4]§. For the sake of concreteness we took Sx = 0.2 and 
Sy = 0.4, which leads to the following frequencies of adatom vibrations: = 1.65 and 
Uy = 2.02. The experimental value for the adsorption energy of K on W is Ed = 2.54 
For the vibration frequency normal to the surface we took = \{uJx + ^y) = 1.84, which 
gives 7 = 0.813. For the interatomic potential (0), we took the parameters Vq = 10 eV 
and Pq = 0.85 A~^. These choices give reasonable values for adsystems0: the interaction 
energies between two adatoms occupying the nearest wells along the furrow and across 
are equal to Vint{asx) ~ 0.98 eV and Vint{cisy) ~ 0.22 eV respectively. Finally, we have 
to define the rate of energy exchange between the adatoms and substrate: we took the 
typical valued r] = 0.1 Ux = 0.165. Note that although some of the parameters are chosen 
rather arbitrary, they are typical for metal atoms adsorbed on metal substrate^i. However, 
as the model is still oversimplified to describe a real adsystem, we have to say that our 
choice of parameters does not claim to provide a quantitative interpretation of the K- 
W(112) adsystem. We do nevertheless believe on the qualitative description of the effects 
under investigation and claim that typical adsystems on anisotropic surfaces (e.g. lithium 
and strontium on molybdenum (112) surface, for which experimental data on the detailed 
coverage dependencies of diffusion characteristics are availableil) should exhibit a similar 



behavior. Finally, for numerical solution of the Langevin equations @, we use the standard 
fourth-order Runge-Kutta method with the time step At = to/20 = 0.19, and the cutoff 
radius was taken as r* = 2asy = 8.94 A. 



III. KINKS 

As the kinks are the main objects of the phenomenological approach!, let us first cal- 
culate their parameters. We recall that kinks can be defined for any commensurate atomic 
structure 60 = s/q, where s and q are relative prime integers; the kink (resp. antikink) 
describes then the minimally possible topologically stable compression (resp. expansion) of 
the commensurate structure. The kink is a quasiparticle, characterized by an effective mass 
nik, a rest energy Sk, and the Peierls-Nabarro (PN) amphtude Epn, corresponding to the 
barrier for the kink translation along the chain. These parameters are determined by the 
dimensionless elastic constant gejf defined as 

9eff = ^Jl ^/nt(«A)- (11) 

Analytically, the kink parameters may be found in the low-coupling limit Qeff ^ 1 or, in 
the strong-coupling (sine-Gordon) limiti geff ^ 1, however, usual real physical systems are 
characterized by the elastic constant Qeff ~ 1, so that both approximations are too crude to 
be applied to our case. For our choice of model parameters, we have g^ff ~ 0.6 for 6q = 1. 
Therefore, we will calculate the kink parameters numerically. 

The numerical method was described in detail in previous paperSil. Briefly, we have 
to choose first an appropriate size of the finite chain in order to insert a single kink into 
the 60 = s/q commensurate background structure; the integers and M must satisfy the 
equation0'0 

qN = sM + a , (12) 

where the kink topological charge a is equal to a = +1 for the kink and a = — 1 for the 
antikink. In the simulation, we restrict ourselves to the concentration range [0.5, 1] because 
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for lower concentrations, the interatomic interaction is too weak and its effects would be 
hardly observable, while at higher concentrations the atoms begin to escape from the first 
adlayer0@. As background structures, we chose the following coverages: = 1/2, = 3/5, 
6*0 = 2/3, 6*0 = 3/4, = 4/5 and Oq = 1. The corresponding values for the number of atoms 
A^o the number of minima of the substrate potential Mq for the every 6*0 are summarized 
in Table |. 

We start with an appropriate initial configuration and allow the atoms to reach the 



minimum-energy configuration (see details in Refs. p6| , |27|) . This determines the structure 



of a kink in its minimal energy state. Then, in order to calculate the parameters that 



characterize the kink translation, we choose a given atom in the kink region (see Ref. |T5D 
and move it to the right by small steps by imposing its x coordinate while all other degrees 
of freedom of the chain remain free to adjust to every new position of the constrained atom. 
This process allows us to find the saddle configuration and therefore, the amplitude of the PN 
barrier Epn as the difference of the saddle and GS energies. Besides, the energy of creation of 
the kink-antikink pair is determined as Epair = -E{kink[0o]} + -E{antikink[^o]} ~ 2£'{GS[^^o]}) 
where E{.} is the energy of the corresponding configuration (notice that it must satisfy the 
relation A^{kink[6'o]} + iV{antikink[6'o]} = 2N{GS[6o]} which is imposed by the total length 
of the atomic chains along x). The kink parameters obtained by this method are summarized 
in Table |I|, and the dependence of the PN energy on the atomic concentration is presented 
in Fig. ||. Note that the function SpniO) has the shape of a devil's staircaseEl. 

From the kink parameters, the phenomenological approach! describes approximately the 
low-temperature behavior of the system as follows. For = 1 and T <^ £pair{Go} Iks the 
concentration of thermally nucleated kink-antikink pairs is equal tc&B 

(^pair) ~ C'exp(-epair/2/cijT) , (13) 

where C ~ (2mfci^^aL/^^B^)"'^^^ s-iid fhk = {nikm^y^'^. For lower coverages 9 = 60 = s/q 
Eq. ( [TB| ) should be properly renormalized, which results in additional factor 1/g in its right- 
hand side. 
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When the concentration 6 shghtly deviates from the commensurate value ^o, the thermal 
kinks are supplemented by the geometrical kinks (if 9 > 9q) or antikinks (if 9 < 9q) with a 
concentration 9geom = (l\9 — 9o\. In the close vicinity of 9q, the total kink concentration can 
be found as 

{9tot) ~ 9geom + '^{9pair)- (14) 

The dimensionless susceptibility x 

X = {0tot)/q^9 (15) 

can then easily be obtained from Eqs. (p!3|-[l^. 

Let us now examine the phenomenology of consequently melted kink superlattice^l. In 
order to describe the atomic mobility in terms of collective excitations, we must first define 
the type of excitations that have to be considered. As an example, let us select a concentra- 
tion 9 in the neighboring of 6*0 = 2/3. For low T, T < epair{2/3}/2kB, we have to use the 
superkinks defined on the background of the 9q = 2/3 structure in the expressions (|T^-|T3|). 
However, in the intermediate temperature range i.e., epair{2/3}/2kB <T < epair{i/2}/2kB, 
when the superkinks are destroyed by thermal fluctuations while the trivial kinks (defined 
on the 6*0 = 1/2 background) are not yet destroyed, we have to substitute the parameters of 
the trivial kinks in Eqs. (p!^- [I^) . In particular, we should take g = 3 for low T but q = 2 for 
intermediate temperatures. In the latter case, however, the parameters of trivial kinks may 
seriously differ from those calculated for the ideal case, because the concentration of trivial 
kinks at the 9 = 2/3 coverage is very large and their interaction is not negligible. 

When the amplitude of the activation barrier for the kink motion is known, the diffusion 
coefficient for a single-kink random walk may be approximately calculated with the 
Kramers theory. For T < Epn/ks this approach gives 

Dk = Dkoexp{-epn/kBT) , (16) 

where 
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if r]l<7] < UJpn, 

if r]> Upn . 



(17) 



Here Upn ^ {■K/asx)i2epn/my'^, a = q- a^x, and r]i = ujpnksT /2'Kepn. 

If the interatomic interaction is strong enough, the inequahty Spn < Spair may easily be 
fulfilled. In this case, within the temperature interval Epn < ksT < Spair-, the kink diffusion 



coefficient is approximately equal to (e.g. see Refs. |33|-p5|) 



nikrj 



1 



I f £ ^2 



8 \knT 



(18) 



Knowing the kink diffusion coefficient, we can find the chemical diffusion coefficient with 
the phenomenological approach! by the formula 



'tot 



(19) 



where for 9 > 9o we should take {9k) = 9geom + impair) and (6'^) = {9pair), while in the 9 < 9q 
case, we have to substitute {9^) = {9pair) and (6*^) = 9geom + impair)- Finally, the chain 
mobility may be found as -B = x^c- These predictions should be now compared with the 
results of simulation and it is the subject of the following section. 



IV. MOBILITY 

To study the mobility, we use an algorithm where we look first for the minimum-energy 
configuration of the system. Then, we increase the temperature up to a given temperature 
T by small steps AT = T/50 during the time ttherm = 100 to = 381. At that point, we apply 
a small dc force F = 0.01 which is gradually increased from F = 0toF = 0.01 during the 
time 100 to, and wait during ty^ait = 100 to in order to allow the system to reach a stationary 
state. Then, for the discrete times tn = n to, we measure the average velocity (v^) of the 
atoms during trun = 100 to, and finally, the procedure is repeated Uave times {uave = 5 in 
the simulation) with different initializations of the random number generator in order to 
estimate the error bars. 
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To demonstrate the effect of the transverse degrees of freedom on the atomic mobihty, 
we considered three different cases of the generahzed Frenkel-Kontorova model: 

• a purely one-dimensional atomic chain with atomic movement restricted to the x di- 
rection (ID), 

• a quasi-one- dimensional atomic chain with two transverse degrees of freedom y and z 
(we will call it the quasi- ID), 

• a true two-dimensional extension of the FK model (2D). 

Note that the interaction between the atoms always has the general form of Eq. (|^ i.e. 
it has a 3D character in the quasi- ID and 2D cases. 

The quasi-lD case can be easily obtained from the general case by choosing a period of 
one lattice constant in the y direction. Namely, we put My = 1 so that all chains move in 
the same way and chose = 5Nq{6} and = 5Mq{6}, where the integers A^^o and Mq are 
taken from Table |I| for each coverage 6 in order to have 5 kinks or antikinks over the length 
under investigation. 

The results of the simulations are presented in Fig. ^ As expected from the phenomeno- 
logical theory, at low temperatures B{6) does have local minima not only for the trivial 
concentrations 6 = 1/2 and 6 = 1, but also at the commensurate concentrations 6 = 2/3 
and 6 = 3/4. These two minima, which involve a kink lattice, disappear when the temper- 
ature is increased, while the minima for the trivial structures survive at any temperature. 
In the simulation, minima do not appear for the other complicated GS structures (e.g. 
6 = 3/5 and 6 = 4/5) because these higher order structures correspond to too low "melting" 
temperatures of the kink lattice. 

In order to check completely the phenomenological theoryi, it would be interesting to see 
if these extra minima appear at very low temperature, but the mobility is then too small to 
allow us to obtain accurate results in a numerical simulation. As the "melting" temperature 
is determined by the magnitude of the effective elastic constant of the kink lattice, one could 
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attempt to increase the parameter Vq in Eq. (^. But in that case, the repulsion between 
the atoms is too large and they begin to escape from the minima of the substrate potential 
in the direction orthogonal to the chaini^. To prevent this escaping and allow the study of 
higher order minima, we artificially restricted the atomic displacements to the x dimension. 
This allowed us to take Vq = 100 eV. This case corresponds to the so-called ID case. The 
results are shown on Fig. |^. 

As seen from Fig. ^, in the ID case the function B{d) has pronounced local minima at 
6 = 2/3 and ^ = 3/4 at much higher temperatures than in the quasi-lD case. One can note 
also that the minimum at = 3/4 disappears when the temperature is increased, while the 
minimum at = 2/3 survives in the whole range of investigated temperatures, since it has a 
greater melting temperature. However, the minima at 6 = 3/5 and 6 = 4/5, where the kink 
lattice has the period 5asx, are not found even for this very strong interatomic repulsion. 
Note also that in real physical systems, such as adsorbed layers, the observation of local 
minima for these complicated structures is unlikely due to existence of transverse degrees of 
freedom. 

In a one- dimensional model at high enough temperature (T > epair/ks), the mobility 
can be calculated with a perturbative approach starting from a system of noninteracting 
atoms. The function B{6) is given 



where Bf = l/mr], qa = do-sx is the average interatomic distance and the elastic constant 
qa is defined by qa = a1^Vll^^{aA) /'^Ti'^esx- Note that the function ( pO]) has local minima 
for trivial configurations only, in agreement with the phenomenological theory in the high 
temperature range. The chemical diffusivity D^. could be obtained replacing the prefactor Bf 
in Eq. ( PUD by a\Vll^f.{a a) / mrj . Figure^ shows that Eq. describes the high-temperature 
simulation results of the purely onedimensional FK model with a good accuracy (except 
in the vicinity of the coverage 6 = 2/3, where epair/ks > T even for the highest studied 




(20) 
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temperature). For this model, the highest possible mobility Bf = l/mr] ^ 6.0, which 
corresponds to the case of noninteracting free atoms, is reached in the middle of the interval 
0.5 <9 < 1.0. 

The high temperature mobility for the quasi- ID chain with transverse degrees of freedom 
(curve (4) in Fig. ^ is approximately two times lower than the values calculated with 
Eq. (pO|) for corresponding parameters. Note that, for the chosen set of parameters in a 
quasi-lD model, T = 0.05 eV is the highest temperature for which the determination of the 
mobility in the first monolayer of atoms is possible. At higher temperatures the atoms start 
to escape from the first layer to the second one, which may seriously distort the results (the 
curve (4) in Fig. ^is not plotted at 6 > 0.8 for this reason). Even if one takes into account 
the fact that the high temperature range required for the validity of Eq. ( pOD is not reached, 
the disagreement between Eq. ( pO]) and the simulation results is large in the quasi- ID case. 
This shows that the presence of the transverse degrees of freedom has the same effect on 
the mobility than an additional friction in the system. This can be understood because 
some part of the work done by the external force is used to excite the transverse degrees of 
freedom. 

Finally, we also simulated the 2D Frenkel-Kontorova system with My = 30, = 
Mo{9} and = 30Nq{9}. The results, presented in Fig. |[ show that there is no essential 
difference between the B{6) dependencies for the quasi- ID and 2D systems except that 2D 
dependencies are systematically lower. The role of the transverse degrees of freedom, already 
noticed for the quasi-lD model, show up again here. It is interesting to notice that Fig. ^ 
shows for the 2D model at T < 0.01 eV, the additional small minimum of B{6) a.t 6 = 4/5 
predicted by the phenomenological theory, which reflects the existence of the kinks/ ant ikinks 
on the background of this coverage. 

The plots of the mobility B versus inverse temperature shown in Fig. |^ for the 2D 
model at selected coverages show that the atomic mobility has an activated character in 
the investigated range of temperatures and coverages. The same qualitative behavior was 
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found for the ID case. Using an Arrhenius form B{T) = Boexp{—Ea/T), we can calculate 
the activation energy Ea and the prefactor Bq. Their dependence versus coverage is shown 
in Fig. ^. The activation energy Ea has a sharp maximum at the coverage 6 = 2/3 which 
corresponds to a well-defined commensurate structure, while activation barriers on both 
sides of = 2/3 are much lower due to the presence of residual kinks/antikinks; the barrier 
for "kink" coverage 6 = 2/3 + 5 is lower than that for "antikink" coverage 6 = 2/3 — 5. 
On the other hand, the maxima of Ea at the higher order commensurabilities 6 = 3/4 and 
6 = 4/5 are much less pronounced. This is consistent with the fact that these higher order 
commensurabilities hardly show up in the mobility curves of Fig. ^. 

The main difference between the quasi- ID system and the true 2D model is due to the 
interactions of kinks in the nearest neighboring channels. For the repulsive interatomic 
interaction studied in the present work, kinks in the nearest neighboring channels repel each 
other at the 6=1 coverage, while for any 6 < 1 they attract each other and tend to form 
domain walls. With the short-range (exponential) interaction studied in our simulations, 
two kinks belonging to neighboring channels attract each other with a potential Vkk{x) oc 
contrary to the usual lawi Vkkix) oc x^. As a consequence, the stiffness of the domain walls 
vanishes and they can be destroyed by thermal fluctuations or external forces for any T ^ 
or F 7^ 0. This is confirmed by the observation of snapshots of the atomic configuration in 
the 2D model. During the time evolution a domain wall of kinks is destroyed as soon as the 
temperature is high enough to provide a noticeable value of the mobility. For example Fig. |^ 
demonstrates the evolution of such a domain wall defined on the background of the 6* = 1/2 
commensurate structure. This case was chosen because its kinks have the simplest structure 
and are more visible than kinks defined on the background of any other more complex 
commensurate structure. The initially well defined kink wall (relaxed configuration for 
T = 0) is smeared out at T = 0.02 eV and F = 0.01, although these values of temperature 
and external force provide a very low value of the atomic mobility {B ^ 0.03) at the 
chosen coverage. This instability of the kink domain walls explains why the true 2D model 
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give results which are not fundamentally different from the results of the quasi- ID case. 
Nevertheless, the interaction between atoms and kinks in the nearest neighbor channels 
does contribute to the dynamics of the system. It results in particular in the lowered values 
of mobility for 2D case in comparison with those for a ID system. However, in more realistic 
2D models with long-range interatomic forces such as elastic or dipole-dipole forces due to 
the substrate, the role of the domain walls might be more essential. 



The chemical diffusion coefficient is more difficult to calculate by MD simulations than 
the mobility. It can be determined in two ways. First, the susceptibility x can be calculated 
with one of the methods described in Ref. ^ and then Df. could be derived from the relation 
Dc = B /x- However this approach relies on the accuracy of the two factors. In the present 
work, we use a direct approach based on the Fick law (^. We start from an nonuniform initial 
concentration profile 6{x) and observe its evolution with time at a given temperature T (we 
will now use the notation 6 instead of ((p)) in the diffusion laws, assuming the existence of 
local equilibrium). The variations of the chemical diffusion coefficient with concentration 
determine the diffusion profile0. For instance, for an approximately constant flux J = 
—Dcdd/dx, fiat sections of the observed concentration profile (low 89 /dx) correspond to 
enhanced diffusivity Dc-, while sharp changes of concentration within some concentration 
interval (high 89 /dx) indicate a lower diffusion coefficient Dc- 

Quantitative data on the variation versus 9 of the diffusion coefficient can be obtained 
by studying the concentration profiles given by the one- dimensional diffusion equation 



The simplest case is the diffusion of an initial step-wise profile in a spatially infinite system 
which gives an explicit expression for the Dc{9) function by the Boltzmann-Matano formula 



V. DIFFUSION 




(21) 



(see e.g. Ref. 0). However, with periodic boundary conditions, computational limitations 
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do not allow us to chose a period large enough to observe such a profile. Therefore, we 
first derive an approximate expression of D^iO) using the phenomenological equations 
|T9|), and then solve Eq. pl| ) with this Dc{6) and periodic boundary conditions. Finally, we 
compare the calculated profiles with those obtained from MD simulation for the same initial 
distribution. 

Let us first apply the kink-gas phenomenology to the determination of Dc. We have 
chosen the room temperature T = 0.025 eV (290 K) since it divides the whole investigated 
coverage interval [0.5, 1.0] into two parts which differ by the mechanism of kink diffusion. 
For the coverage range 0.5 < 9 < 0.66 where the condition T < Spn/kB is satisfied (see 
Fig. |l|), the diffusion of kinks has an activated character and, in terms of the Arrhenius 
representation of the chemical diffusion coefficient D^. = DQexp{—Eac/T), this means that 
the Dc{6) dependence is determined mainly by the variations of the activation energy (note 
that Eac ~ £pn according to the Eq. (|1^)). On the other hand, for 0.66 < 9 < 1.0, where 
T > Spn/kB and where the Peierls-Nabarro energy Spn shows only minor changes with 
coverage, the mass transport is carried out by free diffusion of kinks and the main variations 
of chemical diffusion coefficient will arise from the prefactor Dq. 



One should keep in mind that the phenomenological equations (|l3| -p!9|) can be used 
only for those kinks which are well defined as quasiparticles for a given temperature and a 
given coverage interval, i.e. the condition T <^ Spair/ks must be satisfied. In other words, 
the concentration of thermally excited kink/antikink pairs (|13D for a given structure = 
s/q cannot exceed the maximal possible value 1/q. For the present study, it means that 
the quasiparticles which should be taken into account at T = 0.025 eV are the trivial 
kinks/antikinks of the trivial GS 9o = 1/2 and = 1 and superkinks defined on the 
background of 9q = 2/3 structure. We pointed out in Sec. |T|that, for our model parameters, 
accurate kink parameters can only be determined numerically. But, in order to solve Eq. ( pT]) 
we need some expression for the kink masses. Since the maximum value of the dimensionless 
elastic constant geff defined by Eq. (|TTp approximately equals 0.6 for the chosen set of model 
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parameters, the best estimate is given by the low-couphng hmit expression!! ^ 
m/q^. 



Theoretically, the application of the Eqs. (|T3|-|19]) for the determination of Dc is only 
strictly valid in the close vicinity of the commensurate structures, where the concentration 
of residual kinks is low (in our case, near the coverages 9q = 1/2^ 9q = 2/?) and 9q = 1). Since 
we need Dc{6) for all intermediate 6 values, we have to interpolate between these specific 6 
values. We calculated the values of Dc (indicated by the plus signs in Fig. ||) up to the middle 
points between these specific coverages and used a weighting coefficient, plotted in inset in 
Fig. to mark the significance of each point in the subsequent interpolation procedure. The 
final form of Dc{6) is taken as a superposition of tanh(^) functions chosen to provide a good 
fit of the value deduced from the phenomenological theory around the coverages 6o = 1/2, 
6o = 2/3 and = 1 where it is accurate. Although this procedure cannot avoid some 
arbitrariness, we keep it to a minimum by putting the weighting factor to zero wherever 
the theoretical formula for Dc{6) is not valid. The general shape of the interpolated Dc{6) 
refiect the expected general variations of the chemical diffusivity versus coverage. We see 
that this function monotonically increases in the region 1/2 < 6 < 2/3 (corresponding to 
the decrease of the activation energy Eac which can be deduced from Fig. 0), while at higher 
coverages Dc start to decrease. The high-temperature behavior of Dc at high coverages is 
given in the kink-gas approach by the variation of kink mass. It can be also interpreted in 
terms of the Arrhenius formula Dc = DQexp{—Eac/T). Generally, Dq and Eac change in 
a similar manner (so-called compensation effecM). At high enough temperatures, the slow 
decrease of the Eac at 9 > 2/3 (see Fig. [|) leads only to a slow change of the exponential 
term of the Arrhenius formula. The fast drop of Dc must thus be attributed to the prefactor 
Do. 



Once Dc is known, the second step is to solve Eq. ( pT]) with this Dc dependence and com- 
pare with the MD simulations. To deduce a local coverage from the MD atomic configuration 
at a given time t, we calculate the occupation numbers n{ix, iy] t) (where = 1, . . . , and 
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iy = 1, . . . ,My) defined as the number of atoms in tlie given elementary cell (ixjiy)- The 
results of the simulations for the 2D model (M^. = 84, My = 60 and = 3780) are presented 
in Fig. ^ We started with an artificially prepared step-wise initial configuration: concen- 
tration ^ = 1 in the central region of the lattice (for = Mx/^ + 1? • • • ? 3Mj;/4) and 9 = 0.5 
outside of this region. The system is then allowed to evolve according to the Langevin equa- 
tions (H). The concentration profile 9{ix,tn) = Y!il=i''^i.'^x,'i'y]tn) / My is recorded at dates 
tn = nto- The simulation was repeated five times in order to average the profiles and de- 
crease statistical fluctuations. The simulation profiles for different times are shown in Fig. ^ 
with symbols. They have a flatter section in the middle of the studied coverage interval, 
corresponding to the coverage region with enhanced diffusivity. Moreover, in the same fig- 
ure, we plot with solid lines the theoretical profiles obtained from the numerical solution of 
the diffusion equation (^) with Dc{9) plotted in Fig. |[ The data presented in Fig. |^ show 
that the theoretical and simulation results are in a very good agreement which validates the 
phenomenological approach used for determining the diffusion coefficient. 



VI. DISCUSSION IN RELATION WITH EXPERIMENTS 

It is important to examine the applicability of the theoretical results to real physical 
systems such as atomic layers adsorbed on crystal surfaces. Although experiments can- 
not provide results as detailed as the numerical simulations, a comparison is possible. Our 
model is oversimplified to describe quantitatively a real adsystem — although we chose some 
model parameters close to those available for the K-W(112) adsystem — , mainly because 
the interatomic interaction in real adsystems is much more complicated than the exponen- 
tial interaction used in the present worlS. In the case of adsorption on isotropic surfaces, 
diffusion is affected by the formation of domain walls, especially when the interatomic in- 
teraction is long-range. Our results are more suitable to describe highly anisotropic surfaces 
for which the interaction between neighboring channels is sufficiently weak to reduce the 
role of two-dimensional domain walls. We obtain a qualitative agreement with experiments 
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on diffusion of atoms adsorbed on highly-anisotropic furrowed surfaces. 

There are very few experimental data on the variation versus coverage of the diffusion co- 
efficient for atoms adsorbed on highly anisotropic (furrowed) surfaces. Some data have how- 
ever been obtained using the field emission fluctuations method^ for the K-W(112) and the 
diffusivity was found to increase strongly in the region of the commensurate-incommensurate 
transition: this was interpreted in terms of fast diffusion of soliton^. 

Detailed dependencies Dc{0) and EadO) in the wide coverage interval [0.05, 1.5] are avail- 
ableii for the Li-Mo(112), where the interaction between Li adatoms on Mo(112) is long- 
range and anisotropic. Besides the short range forces, the interaction between the adatoms 
includes also a dipole-dipole repulsion and an oscillating part due to substrate-mediated 
electron exchanged. This is responsible for the existence of peculiar chain-like structures 
p(lx4) and p(lx2), formed by first-order transitions, for coverages In this range 

of 9, the diffusivity Dc was found to depend only weakly on coverage. At higher coverages 
9 > 0.5, the repulsion between Li adatoms starts to play a larger role. This results first 
in the formation of one-dimensional incoherent structures at ~ 2/3, and then the adlayer 
exhibits a one-dimensional compression along the direction of furrows0. In this coverage 
range, Dc{9) was found to increase strongly and monotonically with coverage at low tem- 
peratures. The sharpest increases of diffusivity at low temperatures (T < 250K) appear for 
the commensurate coverages 9 = 2/3 and 9 = 1. This behavior of Dc coincides qualitatively 
with the predictions of the kink-gas approach! and our numerical simulations. Moreover the 
activation energy Eac for chemical diffusion, obtained in Ref. ^ from the slopes of Arrhenius 
plots of Dc, exhibits a monotonical decrease as coverage increases (except for small maxima 
at coverages slightly above the commensurate values 9 = 2/3 and 9 = 1), which is close 
to the behavior of Spn in Fig. |I|. It is also interesting to notice that if Dc{9) is plotted at 
temperatures higher than 300K from the values of prefactor the Dq and activation energy 
Eac measured in Ref. ^ it shows a nonmonotonic behavior with a minimum around 9 = 1, 
which is very similar to the behavior of Dc a.t T = SOOi^ in the two-dimensional FK model 
considered here (see section 0) . 
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Finally, preliminary results of the diffusion study in Sr-Mo(112) systema have demon- 
strated, that diffusivity increases sharply at coverage 9 0.5 which correspond to the com- 
mensurate (4x2) structure of strontium atoms (while at higher coverages Sr atoms form 
incommensurate structures) . One may speculate that this enhanced diffusivity is provided 
by the fast kink diffusion. 



VII. CONCLUSION 

The aim of this paper was to check the predictions of the kink-gas approach with the 
help of Molecular Dynamics simulations. The validity of the kink gas approach ought to be 
questioned because it is based on one-dimensional models while real adsystems are 2D (or 
even 3D taking into account the possible motion of adatoms orthogonal to the surface). 

First, we compared the kink-gas approach [Eqs. (|T^-[1^)] and high-temperature formula 
(pop for mobility B and chemical diffusivity D^. = ksTB/x with the results of our simulation 
of FK models with transverse degrees of freedom. We have found only a qualitative agree- 
ment between the B{6) dependencies obtained in MD simulation of 2D and ID model with 



transverse degrees of freedom and those predicted by kink-gas approach. In section |^ we 
showed that the mobility B is strongly reduced when additional transverse dimensions are 
involved in the system. A quantitative estimation of B using Eqs. (p!^-|T9D shows that the 
kink-gas approach overestimates significantly the mobility unless we artificially introduce a 
higher effective friction ?7e//- This can be understood because some part of energy brought 
into the system by the external force is absorbed by extra degrees of freedom. For example, 
in the simplest case of harmonically interacting atoms at high temperatures, the mobility 
has to be renormalized by a factor 1/3 (i.e. rjeff = 3r]) due to the presence of 2 additional 
degrees of freedom, since the energy is redistributed uniformly (~ ksT /2) between all three 
degrees of freedom. But in our case, where interactions between atoms is anharmonic and 
mobility is investigated at low temperatures, a reliable determination of effective friction 
Tjeff is not possible. 
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By contrast, the study of chemical diffusion (Sec. 0) demonstrates qualitative and quan- 
titative agreement between MD simulation data and the kink-gas approach. The reason is 
that in the case of thermal diffusion, the chemical diffusion coefficient Dc = ksTB/x does 
not change with additional degrees of freedom since the "external force" is provided by the 
thermal energy of the system and is proportional to the gradient of chemical potential /i. 
As the Fick law @) may be rewritten as J = pBV^ = pB{dfi/dp)V p = kBTB/xVp, where 
X = {d In p)/ {dp/ kBT), it is clear that an additional transverse dimension to the system 
leads simultaneously to an increase of free energy and chemical potential p. In other words, 
the decrease of system's mobility due to an extra transverse dimension of the system is com- 
pensated by a corresponding decrease of susceptibility x of the system, so that Dc remains 
approximately the same. 

Our results allow us to conclude that the phenomenological kink-gas approach provides 
a good qualitative explanation not only for Molecular Dynamics simulation data of the mo- 
bility and diffusivity versus atomic coverage in the generalized 2D Frenkel Kontorova model, 
but also for some experimental results on the coverage dependence of surface diffusion. Ob- 
viously, for a better description of the real adsorption systems, the Frenkel Kontorova model 
should take into account long-range, anisotropic, realistic interatomic interaction and the 
presence of surface defects. 
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TABLES 

TABLE I. Parameters of kinks: Nq number of atoms, Mq number of minima of the substrate 

potential for one period of tlie system along x; £pair creation energy of a kink-antikink pair and 
£pn amplitude of the Peierls-Nabarro potential. 

structure Nq Mq Spair (eV) Epn (eV) 

antikink[l/2] 21 43 — 0.378 

6*0 = 1/2 21 42 0.759 — 

kink[l/2] 21 41 — 0.0849 

antikink[3/5] 22 37 — 0.0848 

6*0 = 3/5 21 35 0.007 — 

kink[3/5] 20 33 — 0.0813 

antikink[2/3] 21 32 — 0.0812 

6*0 = 2/3 20 30 0.170 — 

kink[2/3] 21 31 — 0.0192 

antikink[3/4] 20 27 — 0.0184 

^0 = 3/4 21 28 0.055 — 

kink[3/4] 22 29 — 0.0087 

antikink[4/5] 19 24 — 0.0086 

6*0 = 4/5 20 25 0.018 — 

kink[4/5] 21 26 — 0.0071 

antikink[l] 21 22 — 0.0071 
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FIGURES 

FIG. 1. Peierls-Nabarro energy versus coverage. Coverages corresponding to simplest commen- 
surate structures are shown with dashed hnes. 

FIG. 2. The mobihty B of the quasi-lD FK model with transverse degrees of freedom as a 
function of the coverage 6 at selected temperatures T = 0.0025 eV [curve (1)], T = 0.005 eV (2), 
T = 0.020 eV (3), T = 0.050 eV (4). 

FIG. 3. The mobility B versus coverage 6 for purely ID model with Vq = 100 eV at different 
temperatures T = 0.005 eV (diamonds and dotted line), T = 0.05 eV (asterisks and dashed line), 
T = 0.10 eV (triangles and solid line). For a better presentation, the data for the two lower 
temperatures are plotted only within 0.72 < 6 < 0.8, because at other coverages they are the same 
as for T = 0.10 eV. The dash-triple-dotted line is the Eq. (M) for T = 0.10 eV. 



FIG. 4. The mobility B versus coverage 9 for 2D model at selected temperatures T = 0.005 eV 
[curve (1)], T = 0.010 eV (2), T = 0.020 eV (3), T = 0.030 eV (4), T = 0.050 eV (5). 

FIG. 5. The mobility B of the two-dimensional model versus the inverse temperature: (a) for 
coverages 9 = 0.51 (diamonds) and 9 = 0.60 (triangles), (b) kinks (squares), antikinks (diamonds) 
and the background commensurate structure (triangles) for = 2/3. (c) for coverages 9 = 0.72 
(diamonds), 9 = 0.80 (triangles). 

FIG. 6. The activation energy Ea and prefactor Bq for the mobility versus coverage 9 in the 
case of the two-dimensional FK model. 

FIG. 7. Snapshot pictures for two-dimensional model at 9 = 0.51. (a) The initial relaxed con- 
figuration (T = 0), (b) The atomic configuration after the evolution time t ~ 1600 at T = 0.02 eV 

FIG. 8. Chemical diffusion coefficient Dc versus coverage for T = 0.025 eV. Crosses correspond 
to the Dc values calculated with phenomenological equations (P^19); the full curve corresponds 



to the Dc{9) dependence, interpolated with the help of weighting coefficient presented in the inset. 
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FIG. 9. Evolution of the coverage profile versus time: t = (a), t = 190 (b), t = 763 (c) and 
t = 1715 (d). The triangles correspond to the simulation of the 2D model with Mx = 84, My = 60 
and N = 3780 at room temperature T = 0.025 eV, and the full curves are the solution of the 
diffusion equation (pi]). Coordinate x is indicated in lattice units asx 



29 



